Improve efficiency of cut-cell MeshData
#167
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## dev #167 +/- ##
==========================================
- Coverage 97.12% 96.34% -0.78%
==========================================
Files 25 26 +1
Lines 3167 3724 +557
==========================================
+ Hits 3076 3588 +512
- Misses 91 136 +45 ☔ View full report in Codecov by Sentry. |
Merged
|
A script to check accuracy of the resulting hybridized SBP operators: using StartUpDG
using LinearAlgebra
cells_per_dimension = 2
circle = PresetGeometries.Circle(R=0.6, x0=0.0, y0=0.0)
N = 3
rd = RefElemData(Quad(), N) #, quad_rule_face=gauss_quad(0,0,N))
@time md = MeshData(rd, (circle, ), cells_per_dimension;
quad_rule_boundary = gauss_quad(0, 0, N+1),
precompute_operators=true)
@profview_allocs MeshData(rd, (circle, ), cells_per_dimension; precompute_operators=false)
@profview MeshData(rd, (circle, ), cells_per_dimension; precompute_operators=false)
# target_boundary_degree = 2 * N^2 + (N - 1)
# rq_boundary, wq_boundary = gauss_quad(0, 0, ceil(Int, (N^2+(N-1)-1)/2))
# Vtarget = vandermonde(Line(), target_boundary_degree, rq_boundary)
# e, f = 1, 1
# (; cut_face_node_indices_by_elem_by_face) = mt.cut_cell_data
# ids = cut_face_node_indices_by_elem_by_face[e][f]
# Vtarget_nxy = [Diagonal(md.nxJ.cut[ids]) * Vtarget Diagonal(md.nyJ.cut[ids]) * Vtarget]
# w_pruned, inds = StartUpDG.caratheodory_pruning_qr(Vtarget_nxy, wq_boundary)
mt = md.mesh_type
(; volume_interpolation_matrices,
face_interpolation_matrices,
differentiation_matrices,
mass_matrices) = mt.cut_cell_operators
(; wJf) = mt.cut_cell_data
wf = wJf ./ md.Jf
(; x, nxJ, nyJ) = md
(; cut_face_node_indices_by_elem_by_face) = mt.cut_cell_data
Qxyh, hybridized_project_interp_matrices,
hybridized_projection_matrices, hybridized_interp_matrices =
hybridized_SBP_operators(md)
# for e in eachindex(Qxyh)
# Dx, Dy = differentiation_matrices[e]
# M = mass_matrices[e]
# Qxh, Qyh = Qxyh[e]
# Vh = hybridized_interp_matrices[e]
# weak_sbp_error = norm(sum(Qxh, dims=2))
# accuracy_error = norm(Vh' * Qxh * Vh - M * Dx) + norm(Vh' * Qyh * Vh - M * Dy)
# @show weak_sbp_error, accuracy_error
# end
(; x, y, xq, yq) = md
u = @. sin(pi * x) * sin(pi * y)
dudx_exact = @. pi * cos(pi * xq) * sin(pi * yq)
dudx = similar(md.xq)
dudx.cartesian .= rd.Vq * ((md.rxJ.cartesian .* (rd.Dr * u.cartesian) +
md.sxJ.cartesian .* (rd.Ds * u.cartesian)) ./ md.J.cartesian)
for e in axes(u.cut, 2)
Dx, Dy = differentiation_matrices[e]
Qxh, Qyh = Qxyh[e]
Vh = hybridized_interp_matrices[e]
Ph = hybridized_projection_matrices[e]
Vq = volume_interpolation_matrices[e]
dudx.cut[:, e] .= Vq * Ph * Qxh * Vh * u.cut[:, e]
end
@show sqrt(sum(md.wJq .* (dudx - dudx_exact).^2))
# scatter(md.xyz..., zcolor = dudx - dudx_exact, leg=false, ratio=1) |
jlchan
added a commit
that referenced
this pull request
May 22, 2024
* Make triangle and tet node orderings consistent (#153) * simplifying test * fixing vertex orderings * splitting wedge-pyr MeshData tests out * fix triangulate tests * fix wedge mesh ordering * updating hardcoded values in VTK tests * update tests * bump compat for NodesAndModes to 1 * Bump Julia and RecursiveArrayTools compat, remove NamedArrayPartition (#159) * remove NamedArrayPartition files d * bump compat * add sparsearrays compat entry tests seem to be failing on CI without it? * bump Julia compat * bump julia CI version * set lower compat bound for SummationByPartsOperators * bump doc Julia version * reexporting NamedArrayPartition * Refactor `RefElemData` (#157) * add comments * comments and slight refactor * comments and renaming Gauss -> TensorProductGaussCollocation * update docstrings * fix comment grammar * remove unnecessary specialization * removing outdated comments * minor reorganization * removing unnecessary functions * reorganizing helper functions * minor comment cleanup * comments for Kubatko SBP * comments * introduce MultidimensionalQuadrature dispatch type * fix hybrid mesh RefElemData * fix tests * simplify SBP RefElemData * update invalidations * improve comments * more descriptive error message for RefElemData constructors combining Tri/Tet/Wedge/Pyr with TensorProductQuadrature. This should be easy to implement in the future (Stroud) * Add new cut cell `MeshData` (#165) * add some temporary test files * fix docstrings * remove cruft * update plot made nicer plot of Caratheodory pruning for Christina's proposal * update file name Makes it mroe clear this file is meant for plotting * fixing cut cell demo for v1.0+ * adding Caratheodory pruning * cleaning up cutcell demo * improve efficiency slightly * refactoring functions * improve comments * more refactoring * generalize map_to_interval * add new version of "generate_sampling_points" * add dispatch to preserve old version of MeshData * remove cruft * add new routines to create a cut cell MeshData with positive weights * refactoring * fix construction of cut cell face node indices previously assumed same number of nodes on all faces * fix precomputation of operators * clean up test of SBP property * add test of SBP property using new MeshData * add MomentFitting dispatch for old cut cell MeshData d * remove outdated todo * changing wJf to wf internally * add face node index array * make face centroid computation more compact * specialize connect_mesh * add new MeshData based on subtriangulations * add Subtriangulation() as a quadrature type d * add quadrature type to CutCellMesh meshtype * test both Subtriangulation and MomentFitting * add some docs * allow specifying the target cut cell quadrature degree via keyword arg * committing scratch testing files before deleting d * removing some scratch test files * add tests for the weak SBP property * fix an error when num_cartesian_cells = 0 * remove scratchpad files * Adding hybridized SBP operators (#166) * move old MomentFitting cut cell code into separate file * fix default target_degree for volume quadrature * formatting and modifying default parameters for cut interpolation nodes * setting default boundary quadrature + renaming variables for clarity * code cleanup * storing cut_face_node_indices_by_elem_per_face * adding hybridized SBP operators on cut cells * adding cut cell hybridized SBP tests * update test tolerances * update hybridized SBP test tolerance * Improve efficiency of cut-cell `MeshData` (#167) * diagm -> Diagonal * use views * preallocate sorting permutation vector p * improve efficiency of caratheodory * remove type instability for PhysicalFrame basis * bump compat for PathIntersections - should be more type stable * fix compat for PathIntersections * add docstring * release <= 3.4 compat bound on RecursiveArrayTools * add some NEWS.md updates * remove SimpleUnpack and usages of @unpack * Remove Requires.jl (#169) * remove Requires statements in SummationByPartsOperators.jl ext * remove requires statements in StartUpDG.jl * add Plots extension * remove unnecessary [extras] section (since we have test/Project.toml) * remove unnecessary compat restriction on SummationByPartsOperators.jl * add TriangulatePlotExt * remove Requires.jl * add to NEWS.md
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Chasing allocations